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ABSTRACT 

We present Sedonu, a new open source, steady-state, special relativistic Monte Carlo (MC) neutrino 
transport code, available at bitbucket.org/srichers/sedonu. The code calculates the energy- and 
angle-dependent neutrino distribution function on fluid backgrounds of any number of spatial dimen¬ 
sions, calculates the rates of change of fluid internal energy and electron fraction, and solves for the 
equilibrium fluid temperature and electron fraction. We apply this method to snapshots from two- 
dimensional simulations of accretion disks left behind by binary neutron star mergers, varying the 
input physics and comparing to the results obtained with a leakage scheme for the case of a central 
black hole and a central hyper massive neutron star. Neutrinos are guided away from the densest 
regions of the disk and escape preferentially around 45 degrees from the equatorial plane. Neutrino 
heating is strengthened by MC transport a few scale heights above the disk midplane near the inner¬ 
most stable circular orbit, potentially leading to a stronger neutrino-driven wind. Neutrino cooling 
in the dense midplane of the disk is stronger when using MC transport, leading to a globally higher 
cooling rate by a factor of a few and a larger leptonization rate by an order of magnitude. We calculate 
neutrino pair annihilation rates and estimate that an energy of 2.8 x 10^^ erg is deposited within 45° 
of the symmetry axis over 300 ms when a central BH is present. Similarly, 1.9 x 10^^ erg is deposited 
over 3 s when an HMNS sits at the center, but neither estimate is likely to be sufficient to drive a 
GRB jet. 
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1. INTRODUCTION 

Neutron star-neutron star (NS-NS) and neutron star- 
black hole (NS-BH) mergers are prime candidates for ex¬ 
plaining observed short gamma ray bursts (sGRBs) and 
their afterglows (see, e.g., Berger 2013 for a recent re¬ 
view). The large amount of extremely neutron-rich mat¬ 
ter and available energy make these systems potentially 
capable of ejecting matter up to ^4 ^ 200 through the 
r-process (e.g., Lattimer & Schramm 1974; Freiburghaus 
et al. 1999; Korobkin et al. 2012; Goriely et al. 2013). 
The thermal and radioactive glow of this ejecta is thought 
to cause largely isotropic (depeding on the distribution of 
dynamical ejecta), observable infrared/optical emission 
lasting on the order of hours to days (Li & Paczyhski 
1998; Metzger et al. 2010; Berger et al. 2013; Barnes 
& Kasen 2013). Observation of this so-called kilonova 
would provide key information about the merger to com¬ 
plement gravitational wave observations (e.g., Metzger 
& Berger 2012; Nissanke et al. 2013; Piran et al. 2013). 
Some observational evidence of such a kilonova has al- 
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ready been suggested for GRB130603B (Tanvir et al. 
2013; Berger et al. 2013) and GRB060614 (Yang et al. 
2015). 

Realistic simulations of merging compact objects need 
to account for general relativity, a hot nuclear equation 
of state, magnetohydrodynamics (MHD), nuclear reac¬ 
tions, spectral and angle-dependent neutrino transport, 
and possibly neutrino quantum effects (e.g. flavor os¬ 
cillations). Neutrinos in particular play an important 
role in determining the dynamics, brightness, and color 
of predicted ejecta emission. Neutrino emission and ab¬ 
sorption modify the electron fraction and specific entropy 
of the material, which in turn determine which elements 
form from the cooling ejecta (Roberts et al. 2011; Ko¬ 
robkin et al. 2012; Wanajo et al. 2014) and the result¬ 
ing photon opacities (Barnes & Kasen 2013; Kasen et al. 
2013). Neutrino irradiation can also drive a thermal out¬ 
flow, generally increasing the amount and electron frac¬ 
tion of ejecta (McLaughlin & Surman 2005; Surman et al. 
2006; Wanajo & Janka 2012; Fernandez & Metzger 2013; 
Just et al. 2015a; Martin et al. 2015; Goriely et al. 2015; 
Foucart et al. 2015), especially in the presence of a cen¬ 
tral hypermassive neutron star (HMNS) (Dessart et al. 
2009; Perego et al. 2014; Sekiguchi et al. 2015; Metzger & 
Fernandez 2014, hereafter MF14). Neutrino-antineutrino 
annihilation may generate large amounts of thermal en¬ 
ergy in baryon-poor regions and remains a possible en¬ 
gine driving the GRB jet (Eichler et al. 1989; Meszaros & 
Rees 1992; Popham et al. 1999; Zalamea & Beloborodov 
2011; Leng & Giannios 2014), though many calculations 
show that the energy production is at best marginally 
capable of powering GRBs (e.g. Setiawan et al. 2006; 
Dessart et al. 2009) 
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Due to the great complexity of this problem, all cur¬ 
rent and past simulation efforts make some level of ap¬ 
proximation or evolve only for very short times to make 
the problem computationally tractable (see, e.g., Faber 
& Rasio 2012; Shibata & Taniguchi 2011 for reviews). 
Neutrinos are ignored altogether in many studies for sim¬ 
plicity or efficiency (e.g., Etienne et al. 2012; Kiuchi et al. 
2014; Bauswein et al. 2014; Bernuzzi et al. 2015; Takami 
et al. 2015). Fernandez & Metzger (2013) approximates 
self-irradiation from the disk as a gray lightbulb arising 
from a ring, with optically thin cooling rates corrected 
for optical depth effects. Various forms of the leakage 
scheme of Ruffert et al. (1996) can be used to more 
accurately treat cooling, heating, and electron fraction 
changes from neutrinos whenever the disk becomes op¬ 
tically thick (Janka et al. 1999; Ruffert & Janka 1999; 
Rosswog & Liebendorfer 2003; Sekiguchi et al. 2011a,b, 
2012; Kiuchi et al. 2012; Deaton et al. 2013; Galeazzi 
et al. 2013; Neilsen et al. 2014; Foucart et al. 2014) al¬ 
though ad-hoc assumptions about the angular distribu¬ 
tion of radiation are still needed to compute neutrino 
absorption (MF14; Perego et al. 2014; Fernandez et al. 
2015a, b). Moving from simple approximations to ac¬ 
tual neutrino transport, Dessart et al. (2009) use New¬ 
tonian multi-group flux-limited diffusion (MGFLD) dur¬ 
ing evolution and multi-group multi-angle transport dur¬ 
ing post-processing analysis. General relativistic two- 
moment neutrino transport with an analytic closure (e.g., 
Thorne 1981; Shibata et al. 2011; Gardall et al. 2013) 
is the state of the art in multi-dimensional simulations 
and currently gives the most accurate approximation to 
full Boltzmann transport of all of the methods employed 
in time-dependent simulations. Initial studies have em¬ 
ployed a gray (energy-integrated) transport scheme with 
general relativity to simulate a merger and remnant (Shi¬ 
bata & Sekiguchi 2012; Sekiguchi et al. 2015; Foucart 
et al. 2015), and Just et al. (2015a) recently simulated an 
axisymmetric remnant disk with Newtonian multi-group 
two-moment transport. 

Unlike the above transport methods, Monte Garlo 
(MG) transport is continuous in space, direction, and 
energy, freeing it from many of the grid effects and distri¬ 
bution function approximations used in other methods. 
Though it is more accurate, it is more computationally 
expensive and has seen more limited use in the past. 
MG transport has been used to study neutrino equilibra¬ 
tion in a static isotropic background (Tubbs 1978) and 
to study transport through static spherically symmet¬ 
ric fluid in the context of core-collapse supernovae (e.g., 
Janka & Hillebrandt 1989; Janka 1991; Keil et al. 2003. 
More recently, Abdikamalov et al. (2012) uses a time- 
dependent MG neutrino transport scheme in static spher¬ 
ically symmetric core-collapse simulation snapshots. MG 
transport has also been used in the context of photon 
transport in la supernova explosions (e.g., Kasen et al. 
2006; Wollaeger et al. 2013; Roth & Kasen 2015) and 
accretion disks (e.g., Ryan et al. 2015). 

In this paper, we investigate the effect of neutrinos 
on the rates of change of the composition and ther¬ 
mal energy of the remnant disk and ejecta using time- 
independent Monte Garlo (MG) neutrino transport cal¬ 
culations. We calculate properties of the neutrino radia¬ 
tion fields to pinpoint the regions of largest error in more 
approximate schemes and proceed to estimate the effect 


this would have on dynamical simulations of compact 
object mergers. We begin by introducing the specifics 
of our neutrino transport scheme and the background 
fluid in Section 2. We proceed to describe the observed 
properties and effects of the neutrino radiation held and a 
comparison to those seen by the leakage scheme in the dy¬ 
namical calculation with a central black hole in Section 3. 
In Section 4, we extend the results to fluid backgrounds 
with a central hypermassive neutron star (HMNS). We 
briefly discuss the effects of neutrino pair annihilation for 
both sets of background data in Section 5. In Section 6, 
we discuss the potential implications of these results on 
the dynamical calculation and the implications for nucle¬ 
osynthesis and kilonovae. In Section 7, we conclude and 
list the main points that can be drawn from our results. 

2. METHODS 
2.1. Background Fluid 

Our background fluid snapshots come from the ax¬ 
isymmetric 2D simulations from MF14 of remnant disks 
modeled after those left behind by a binary neutron star 
merger. The disks have a mass of O.OAMq and circle a 
3Mq BH or HMNS. Table 1 lists the times and global 
properties of the fluid snapshots we use. Figure 1 shows 
the background density, temperature, electron fraction, 
and fluid speed at 3 ms after the start of the dynami¬ 
cal simulations. To avoid large inconsistencies in using 
a relativistic transport code with fluid velocities calcu¬ 
lated with a Newtonian code, we cap the fluid speeds at 
a maximum Lorentz factor of 2. 

Though the fluid moves around on an orbital timescale 
of about 3 ms (at R = 50 km), the time required for the 
disk to significantly change its structure is set by the vis¬ 
cous timescale of about 1 s (Fernandez & Metzger 2013). 
The diffusion time for the neutrinos can be approximated 
by tdiff ^ ^ ~ 10“^s, where r ^ 3 is a representative av¬ 
erage optical depth, I ^ 10^ cm is an approximate charac¬ 
teristic size of the portion of the disk opaque to neutrinos, 
and c is the speed of light. The neutrino diffusion time 
is much shorter than the viscous timescale. Hence, al¬ 
though the fluid orbits faster than neutrinos can escape, 
we can safely assume that the global disk structure is 
effectively static on the neutrino propagation timescale. 

The 3Mq central object is well above the maximum 
mass of current observational and theoretical constraints 
for the maximum mass of a cold neutron star (e.g., Lat- 
timer & Prakash 2000; Steiner et al. 2013). The cen¬ 
tral object in these simulations is likely to become a 
black hole, but the temporary stability depends on ro¬ 
tation, thermal support against collapse (e.g. Kaplan 
et al. 2014), and magnetic held strength and configura¬ 
tion, though the latter requires magnetic fields on the 
order of 10^^ G for the effect to be relevant (e.g., Gardall 
et al. 2001). As did MF14, we remain agnostic to how 
long the mass in the central object is supported against 
collapse and consider both instant black hole creation 
and indefinite HMNS stability to bracket the parameter 
space. 

2.2. Opacities and Emissivities 

Neutrino-fluid interactions depend on the fluid density 
p, temperature T, and electron fraction Ye, as well as the 
neutrino energy and species s^. The interactions are 
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Figure 1. Fluid Backgrounds ai t = Sms from the start of the dynamical simulations of MF14. Left: central BH Right: central 
HMNS. The outer radius on the plot is at 250 km and the inner radius is 30 km. Each quadrant covers half of the simulation 
domain. On the top left of each plot is the density, which peaks at around 5 x 10^® g cm“^. On the top right is the magnitude of 
the velocity in units of c. On the bottom right is the electron fraction, where grey colors indicate electron fractions larger than 0.5 
due to very low-density hydrogen that is present for numerical reasons. On the bottom left is the temperature, which peaks around 
5 MeV. The black/white curve is the p = 10° g cm“^ contour, below which Sedonu opacities and emissivities are set to zero (see Section 2.2). 


Table 1 

Background Global Properties 


Time 

(ms) 

■Wdisk 

(Me) 

(T) 

(MeV) 

{Ye} 

^visc 

(B s-i) 

d^core (B S 
each Ue , ide 

(B s-1) 

<^> 

(s-‘) 

(Me) 

^Si (B 

s-') 

(MeV) 

Central 

0 

BH 

3.01(-2) 

2.95 

0.102 

3.02(1) 


5.19(1) 

2.49(1) 

2.24(-6) 

1.58(1) 

4.01(1) 

17.3 

3 

2.93(-2) 

2.81 

0.121 

2.84(1) 


3.29(1) 

2.47 

2.27^7) 

1.54(1) 

1.92(1) 

15.1 

30 

1.52(-2) 

1.93 

0.114 

7.95 


5.77 

4.92(-l) 

1.35^8) 

2.86 

3.19 

10.8 

300 

4.48(-3) 

6.29(-l) 

0.190 

3.63(-l) 


1.93(-2) 

1.47(-1) 

0.00 

4.93(-2) 

1.44(-2) 

5.08 

Central 

0 

HMNS 

3.01(-2) 

2.95 

0.102 

3.05(1) 

1.73(1) 

3.58(1) 

3.33(1) 

2.42(-3) 

1.60(1) 

3.90(1) 

15.4 

3 

3.01(-2) 

2.93 

0.137 

2.89(1) 

1.73(1) 

2.77(1) 

6.51 

2.22(-3) 

2.67(1) 

2.64(1) 

16.0 

30 

3.01(-2) 

2.75 

0.177 

2.11(1) 

1.02(1) 

1.87(1) 

3.14(-1) 

6.46(-3) 

2.59(1) 

1.84(1) 

16.6 

300 

3.01(-2) 

1.09 

0.264 

5.85 

3.22 

4.75 

2.43(-l) 

5.88^3) 

4.99 

3.84 

13.6 

3000 

1.32(-2) 

4.66(-2) 

0.300 

1.19(-1) 

1.02 

-8.41 (-4) 

1.41(-3) 

9.62^3) 

3.98(-7) 

2.73(-6) 

2.43 


Notes: Quantities extracted from the dynamical simulations conducted by ME 14. The numbers in parentheses indicate the power of 
10 with which the data given must be scaled, e.g., 6.95(—1) is 6.95 x 10“^. Mdisk is the mass remaining in the disk. (T) and (Ye) are 
mass-weighted averages of the disk temperature and electron fraction, respectively. F^visc is the integrated viscous heating rate. Ccore 
is the luminosity of each neutrino species emitted from the core. Each of the previous quantities are taken as input to Sedonu, and the 
following can be compared to the Sedonu results. Cv> — LLv is the net rate of energy loss from the fluid by neutrinos and {dYe / dt) is the 
mass-weighted average of the rate of change of the electron fraction computed by neutrino leakage. Mjy is the mass in which neutrinos 
are a larger source of heat than viscosity is, i.e. —Cv> Hvisc- >Csi is the luminosity of the disk self-irradiation assumed to be emitted 
from two rings above and below the equatorial plane. is the energy density-weighted average energy of these emitted neutrinos. 

1 B = 10^^ erg. 


taken into account via neutrino opacities and emissivities 
calculated by NuLib^ (O’Connor 2014) and are output in 
tabular form for a range of values of {p^T^Ye, . 

We use a table spanning p = g cm“^ with 82 

logarithmically-spaced points, T = 0.05 — 200 MeV with 
65 logarithmically-spaced points, = 0.035 — 0.55 with 

® open source, available at www.nulib.org 

input flies available at bitbucket.org/srichers/sedonu 


51 linearly-spaced points, and Ej, = 0.5 — 200 MeV with 
48 logarithmically-spaced bins. We demonstrate that 
this table has sufficient resolution in Appendix A. The 
opacities and emissivities below the table minima in 
{p^T^Ye^Ey} are very low or affect a very small amount 
of mass. They are hence dynamically unimportant, so 
we assume them to be zero. Heavy lepton neutrinos play 
a relatively minor role (see Sections 3 and 4) since they 
deposit energy only via neutral current reactions, so we 
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Table 2 

Included Physics 


Process 

Leakage 

Full 

Shen 

LS220 

NoWM 

NoScat 

NoPair 

NoRel 

Simple 

z/g/^e Emis/Abs on 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

Weak Magnetism Correction 


/ 

/ 

/ 


/ 

/ 

/ 


Elastic Scatter on e, n, p, a 


/ 

/ 

/ 

/ 


/ 

/ 


Ux Emis/Abs/Pair Production 


/ 

/ 

/ 

/ 

/ 


/ 


Special Relativity 


/ 

/ 

/ 

/ 

/ 

/ 



Elastic Scatter on Heavy Nuclei 



/ 

/ 






Equation of State 

Helmholtz 

Helmholtz 

Shen 

LS220 

Helmholtz 

Helmholtz 

Helmholtz 

Helmholtz 

Helmholtz 


Notes: Included physics in each class of neutrino transport simulations. The first column represents the physics included in the leakage 
scheme in the original dynamical simulations by MFI4. Full represents our most complete set of physics, while Simple is designed to 
replicate the physics used in the leakage calculations as closely as possible. The inclusion of the first four processes and the choice of 
equation of state go into generating the NuLib opacity tables. Special relativistic physics is turned on or off within Sedonu. 


simulate three effective neutrino species Sj, = 
where Ux accounts for 

We experiment with excluding various processes and 
corrections as listed in Table 2 to determine how much 
each approximation affects the resulting neutrino radi¬ 
ation field and fluid source terms. The Full simulations 
embody our most complete set of physics which will serve 
as the standard for comparison. The Simple simulations 
account only for charged-current interactions on free nu¬ 
cleons in order to match as closely as possible the physics 
assumed in the leakage scheme used in the dynamical 
simulations of MF14. The Full simulations include weak 
magnetism and recoil corrections and the opacity for each 
neutrino species to scatter elastically on electrons, neu¬ 
trons, protons, a particles, and heavy nuclei, unlike the 
Simple and Leakage simulations. Opacities for scattering 
on heavy nuclei are corrected for ion-ion correlations, the 
heavy-ion form factor, and electron polarization, though 
heavy nuclei are ignored when using the Helmholtz equa¬ 
tion of state (EOS, see below). Each of these processes 
is implemented as described in Burrows et al. (2006). 

Absorption opacities are converted into emissivities 
via Kirchhoff’s Law and account for final-state electron 
and positron blocking. There are additional approxi¬ 
mate emissivities calculated for pair processes, namely 
electron-positron annihilation (e“ -h ^^1 + ^ 1 ) and 
nucleon-nucleon Bremsstrahlung (ni + 77.2 ris -h 744 + 
Ui + Vi where rij represents any nucleon). While it is 
in principle incorrect to apply Kirchhoff’s Law to these 
emissivities to get opacities since the opacities depend on 
both neutrinos and anti-neutrinos (which are not neces¬ 
sarily in equilibrium), doing so gives correct absorption 
rates in the optically-thick and trivial neutrino-free lim¬ 
its. In this way, NuLib yields annihilation rates that are 
correct to an order of magnitude, though it somewhat 
overestimates the effective opacity. In light of this, we 
include pair processes only for heavy lepton neutrinos 
during the transport step (post-processing annihilation 
calculations are described in Section 2.3.3). 

The opacities are corrected for final-state blocking and 
in general depend on the chemical potentials of the par¬ 
ticles involved in reactions (Burrows et al. 2006). The 
nucleon and lepton chemical potentials at a given den¬ 
sity, temperature, and electron fraction depend on the 
details of the equation of state (EOS). To compare as di¬ 
rectly as possible with the dynamical simulations, we use 
the Helmholtz EOS (Timmes & Swesty 2000) including 
neutrons, protons, and a particles in nuclear statistical 


equilibrium (NSE). We also use two popular hot nuclear 
equations of state: those from Shen et al. (2011) and 
the 220-MeV incompressibility version from Lattimer & 
Swesty (1991)^^. 

The opacities calculated by NuLib are given in the rest 
frame and are transformed into the lab frame according 
to Mihalas & Mihalas (1999), 

^lab = ^rest^^ : (1) 

c 

where 7 is the fluid Lorentz factor, v is the fluid velocity 
vector, and D is the neutrino direction unit vector. To 
test the effect of neglecting relativistic transformations, 
we simply set all fluid velocities to 0 . 

2.3. Monte Carlo Neutrino Transport 

Sedonu^^ is a special-relativistic Monte Carlo (MC) 
neutrino transport code, evolved from the photon trans¬ 
port code Sedona (Kasen et al. 2006). Sedonu sim¬ 
ulates neutrinos passing through and interacting with 
a static background fluid snapshot. We import the 
fluid grid structure, density p, temperature T, and elec¬ 
tron fraction Ye from grid-based simulation snapshots of 
any geometry in zero to three dimensions. In a zero¬ 
dimensional background (i.e., a one-zone model), sym¬ 
metries are imposed in all three directions to create a 
completely isotropic and homogeneous background. 

The neutrinos move around in a separate superimposed 
Cartesian space so that the geometry of the underlying 
fluid grid comes into play only when the neutrinos re¬ 
quire information about the fluid at their current loca¬ 
tion. The neutrinos are discretized into packets, each 
of which is specified by the energy Ejy of each neutrino 
in the packet, the location x of the packet, the direc¬ 
tion unit vector D of travel, and the total energy Ep of 
the packet. When special relativity is included, the grid 
structure is assumed to be in the lab frame and the fluid 
properties are given in the rest frame. 

In a real merger, general relativity would diminish the 
energy of outgoing neutrinos and increase the energy of 
incoming ones. To estimate the magnitude of this effect, 
we can assume a Schwarzschild metric outside of and 
sourced only by the 3Mq central object, which implies 
that 

^ n-2GMIr2c‘^V^^ 

\l-2GM/ric^J ’ ^ ’ 

both available in tabular form at stellarcollapse.org 
open source, available at bitbucket.org/srichers/sedonu 
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where M is the mass of the central object and 
and Ejy ^2 are the energies of a given neutrino at radii 
ri and r 2 , respectively. The strongest redshift effect 
we could expect is the difference between the neutrino 
energy at the inner boundary (30 km) and the outer 
edge of the disk 250 km), which comes out to be 
^zy, 250 km/^zy, 30 km = 0.86. The ucutrino opacities scale 
approximately as ^ (e.g., Burrows et al. 2006), re¬ 

sulting in about a 25% effect on the opacities over this 
distance. However, most of the neutrino energy is emit¬ 
ted and absorbed over distances of tens of kilometers, so 
errors from excluding gravitational redshift will be nec¬ 
essarily smaller than this. In general relativity, neutrinos 
would also follow null geodesics rather than straight lab- 
frame lines, but we defer to other authors to determine 
the importance of this (see Section 6 for a discussion). 

By their very nature, MC simulations output data with 
random fluctuations that decrease with the number of 
MC elements. To keep the fluctuations of global quan¬ 
tities in Table 3 below 0.1% we propagate 2 — 4 x 10^ 
particles in each simulation. 

2.3.1. Emission 

Neutrinos are either emitted from the fluid itself, or 
from an off-grid source like a central hypermassive neu¬ 
tron star (HMNS). If the source of neutrinos is the on- 
grid fluid, the number of neutrino packets spawned in 
each grid cell is proportional to the cell’s total rest-frame 
emissivity. Within each cell and in the rest frame, the 
direction and position distribution of spawned neutrinos 
are isotropically randomly determined. The net comov¬ 
ing luminosity of a grid cell is determined by numerically 
integrating the emissivity over energy bins according to 

£ceii = 47ry E E SiAEjy^i , (3) 

species i 

where V is the grid cell volume, Si is the emissivity (units 
of erg s“^cm“^Hz“^sr“^), and AEy^i is the width of the 
neutrino energy bin i. Each emitted neutrino packet in 
a given grid cell has the same packet energy, determined 
by 

E, = , ( 4 ) 

where At is the length of the time step (arbitrary for 
steady-state calculations like these, as it always cancels 
out). Pemit is the number of neutrino packets to be emit¬ 
ted from the grid cell and is proportional to the cell’s 
neutrino energy emission rate according to 

Pemit = PtotalV^- ^ - 5 ( 6 ) 

^cells^^cell 

where Ptotai is the total number of neutrino packets used 
in the simulation (2 x 10^ in our case). Though each cell 
has a different velocity and hence the lab-frame emissiv- 
ities are modified. Equation 5 is used only to distribute 
computational resources through the disk and has no 
physical meaning. 

Neutrino energies fall into one of 48 energy bins match¬ 
ing the NuLib tables we use. The energy bin of any given 
neutrino is chosen by randomly sampling the local neu¬ 
trino energy-dependent emissivity. Neutrinos are emit¬ 
ted only from the center of a neutrino energy bin Ey^i. 


This is to better maintain consistency required by Kirch- 
hoff’s Law between the emissivity of a grid cell and the 
product of the opacity and the neutrino blackbody func¬ 
tion, both of which are also evaluated at the bin center. 
As neutrinos move through fluid cells with different ve¬ 
locities, they are Lorentz transformed away from the bin 
centers, reducing the level of consistency, but we ignore 
this minor discrepancy. 

Once an MC particle is created, it is Lorentz trans¬ 
formed into the lab frame. Additionally, it should be 
noted that the fluid in a moving cell is length contracted, 
such that its rest-frame volume is larger than its lab- 
frame volume by a factor of the Lorentz factor. Thus, 
including special relativity increases the rest mass (by at 
most 4% in any snapshot), average neutrino energy, and 
the net luminosity of moving grid cells. 

In this paper, the temperature and luminosity of elec¬ 
tron neutrinos and electron anti-neutrinos emitted from 
a central HMNS are taken directly from ME14. There, 
Ty^ = 4MeV and Ty^ = 5MeV, and the luminosity of 
each species obeys 


20 B s 


-1 


(30k) t > 10 ms . 


(6) 


The values of the HMNS luminosity at each of our snap¬ 
shots is also listed in Table 1. When heavy lepton neutri¬ 
nos are included, we choose their temperature and lumi¬ 
nosity to be the same as those of electron anti-neutrinos. 
The neutrinos emitted from a central HMNS are given 
an isotropically random position on the inner boundary 
sphere and an isotropically random direction within the 
outward 27r steradians. Neutrino energies are sampled 
from a Eermi-Dirac blackbody distribution of the ap¬ 
propriate temperature and zero chemical potential. The 
HMNS emits 2 x 10^ packets in addition to the 2 x 10^ 
emitted from fluid in the disk, and the energy of each 
HMNS-emitted neutrino packet is then chosen such that 
the total HMNS luminosity is equal to Ecore- 


2.3.2. Propagation 

Once all neutrinos have been created, their motion is 
in the lab frame along straight lines until they escape 
or are destroyed, punctuated by moments of scattering. 
The distance moved along any straight-line segment is 
the minimum of the following computed distances along 
the packet’s direction of travel: (1) the distance to the 
simulation outer boundary ^boundary, (2) dceii, which is 
0.4 times the length of the smallest dimension of the 
cell currently occupied, and (3) the interaction distance 
dinteract- Wc usc this method of calculating dceii rather 
than computing a geometric distance to the cell bound¬ 
ary for efficiency, but we demonstrate that the factor of 
0.4 is small enough to adequately substitute for a more 
precise geometric calculation in Appendix A. The inter¬ 
action distance is randomly sampled such that its prob¬ 
ability density obeys 

^(rfinteract) = , (7) 

where n = t^s E fhe sum of the scattering and ab¬ 
sorption opacities calculated by NuLib. The particle is 
then moved a distance d = min{dboundary, dinteract, c^ceii} 
and interacts with the fluid if the shortest distance was 

dinteract• 
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When a packet interacts with the fluid, it is randomly 
decided whether it scatters or is absorbed, with probabil¬ 
ities of each proportional to their respective opacities. If 
the packet scatters on the fluid, it is transformed into the 
fluid rest frame, given an isotropically random new di¬ 
rection with the same energy, and transformed back into 
the lab frame. When a packet is absorbed, it is simply 
destroyed. 

Irrespective of which of the three above distances was 
smallest, the neutrino packet deposits thermal energy 
continuously into the grid cells it passes through ac¬ 
cording to ^dep = EpptZad and lepton number according 
to A^dep = evaluated in the comoving frame, 

where / = 1 for I = —1 for Pg? and / = 0 for 
Though these quantities may be more or less than the 
actual packet energy and lepton number, depending on 
whether the packet moves more than or less than one op¬ 
tical depth, energy and lepton number conservation are 
well approximated when averaging over many particles. 
See Appendix B.l for a test of this algorithm. 

The lab-frame rate of change of comoving-frame inter¬ 
nal energy density e and electron fraction Tg in n given 
grid cell is then 


= ^ WAt 1 + T ^dep 

\ steps 

\ steps / 

V is the grid cell’s volume in the comoving frame. At is 
the (arbitrary) emission time interval in the lab frame, 
rUn is the mass of a neutron, and T’emit and A^emit are the 
sum of the emitted comoving-frame neutrino energy and 
lepton number, respectively, from all neutrino species. 
The sum is over all steps (propagation segments between 
emission, scattering, absorption, or escape) for all neu¬ 
trino packets in the cell. The specific heating rate due 
to viscosity in the comoving frame is is taken from 
the simulations of MF14 and the fluid Lorentz factor 7 
transforms the time derivative in the viscous heating rate 
to the lab frame. 

Each grid cell has a distribution function con¬ 
sisting of 6144 energy/direction bins composed 
of (8 latitudinal bins) x (16 longitudinal bins) x 
(48 energy bins). The latitudinal bins have constant 
size in cos(^), where 0 is the angle from the pole, so 
each bin covers the same solid angle. As each neutrino 
propagates it contributes an energy Erad = Epd/cAt 
to its corresponding energy and angular distribution 
function bin {(pi^Oj^Uk} in its current grid cell. The 
energy bins match those of the NuLib table. The total 
neutrino energy density Cy and average energy {Ep) in 
each cell are then 

^>^ = ^EEEE^-d (9) 

(pi Oj Vk steps 

= TV E E E E ^-d. (10) 

^ cpi 9j steps 



2.3.3. Annihilation 


In a separate post-processing step that does not feed 
back into the neutrino distribution, we calculate an anni¬ 
hilation rate (erg s“^cm“^) in each cell following Ruffert 
et al. (1997). We integrate an annihilation kernel over 
the distribution according to 


Qi.„ = EEEE 




Eyh - ^rrieC^ 


EijiE. 


n.k 


^ ^u,ij^iy,klEjl(^Eiy i^ Ely COS 0ji^ , 


( 11 ) 


Here, and are the neutrino and anti-neutrino 
contributions, respectively, to the cell’s neutrino energy 
density (erg cm“^) from the neutrino energy bin center 
Ej^^i and direction bin j. The angle between the centers of 
direction bins j and I is Oji- This must then be summed 
over each of the three neutrino-anti-neutrino pairs to get 
the total annihilation rate. Since we group all four heavy 
anti/neutrino species together in our simulations, we set 

before performing the an¬ 
nihilation calculations. See Appendix C for details on the 
annihilation kernel Rji{uicosO^i)- The derivation of 
annihilation rates assumes that the sum of the neutrino 
energies is much larger than the sum of the rest masses of 
the produced electron and positron. Hence, we subtract 
the mass energy of the electron-positron pair from the 
annihilation rate to under-emphasize energy contributed 
near the minimum-energy limit. Additionally, this causes 
the annihilation rate to represent only the deposited ther¬ 
mal energy without counting mass energy. To check how 
large of an effect this has, we calculate the integrated 
annihilation rate at the 3 ms snapshot for both the BH 
and HMNS cases within 45° of the axis of symmetry 
without subtracting the electron rest mass. This caused 
the energy deposition rate from neutrino annihilation to 
increase by only 2.5%. 


2.3.4. Equilibrium 

After all particles have propagated through the fluid 
and have left a tally of how much energy and lepton 
number they deposited in each cell, we can determine 
what combination of {T, Yg} causes the fluid in each cell 
to emit the same amount of energy and as many leptons 
as it absorbed. The equilibrium T and/or Yg is con¬ 
verged upon using Brent’s method (Brent 1973), which 
queries the NuLib tables for emission rates at successive 
guesses of {T, Yg} until the integrated NuLib emissivities 
(both energy and lepton number) match the absorption 
rates calculated during the MC transport. The equilib¬ 
rium values are physically sensible quantities only where 
the timescales for such an equilibrium to be reached are 
short compared with the dynamical timescale. The pro¬ 
cess of transporting neutrinos and solving for equilibrium 
can be done iteratively, allowing temperature and elec¬ 
tron fraction changes to affect the neutrino sources, un¬ 
til a truly time-independent equilibrium is reached, as is 
done in the irradiation tests in Appendix B. 2 . The true 
time-independent solution of the NS-NS post-merger disk 
problem we study here is trivially a zero temperature 
disk, but to evaluate how strongly the fluid and neutrino 
radiation fields are out of equilibrium we stop after a 
single iteration to arrive at a local rather than global 
equilibrium. 
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2.4. Neutrino Leakage 

For completeness, we review the neutrino leakage 
scheme used by Metzger & Fernandez (2014). Through¬ 
out the following, only electron neutrino and anti¬ 
neutrinos are included. 

The central HMNS emits neutrinos with the same tem¬ 
perature and luminosity as described in Section 2.3.1, 
and the neutrino flux due to the HMNS at any given 
location is attenuated by the (grey) optical depth inte¬ 
grated radially from the HMNS. 

Metzger & Fernandez (2014) determine the rate of en¬ 
ergy loss at any location in the torus by interpolating 
between the optically thin free-streaming limit and the 
optically thick diffusion limit, given by the effective lu¬ 
minosity^^ 

oo 

= E 1^//, ^ttV [. (12) 

t + Miff/tloss J 

* 0 

Here, Si is again the the neutrino emissivity for species i, 
iioss = is the characteristic time for the fluid 

to lose its internal energy via neutrino emission, and 
^diff = X {d/c) is the characteristic time for neutri¬ 
nos diffusing over a characteristic escape distance. The 
first term in the expression for tdiff represents a typical 
optical depth through which neutrinos would need to dif¬ 
fuse to escape, where n is the energy-averaged neutrino 
absorption coefficient due to charged-current reactions 
at the given location (see Fernandez & Metzger 2013 for 
details). The second term is the time required for an 
unimpeded neutrino to cross the same distance. The es¬ 
cape distance is taken to be d = min{r, i7_L, 77||}, where 
and are the vertical and horizontal scaleheights, 
respectively. 

The disk self-irradiation scheme assumes that neutri¬ 
nos are emitted from two rings, one above and one below 
the midplane. The location of the ring is at the effec¬ 
tive luminosity-weighted average radius and polar angle 
in each hemisphere, and the luminosity of each ring is 
half of the volume-integrated effective neutrino luminos¬ 
ity. Neutrinos are emitted from both rings with a zero 
chemical potential blackbody spectrum, the temperature 
of which is the effective luminosity-weighted average fluid 
temperature. The fluxes of each neutrino species at a 
given location are independently attenuated by an op¬ 
tical depth Tirr = max(/i:^d,/^i^ringdring), where dring is d 
evaluated at the ring’s location and /i:i,ring is the absorp¬ 
tion coefficient for species i at the ring’s location. For 
details, see Fernandez & Metzger (2013) and Metzger & 
Fernandez (2014). 

Comparing Equation 3 (used in computing >Cemit 
Tables 3 and 4) with Equation 12 (used in computing 
the self-irradiation luminosity £si in Table 1), it is clear 
that £si < >Cemit- The quantities both represent neutrino 
radiation coming from the disk itself in some capacity, 
but Esi is diminished by optical depth effects, preventing 
direct comparison with £emit- 

3. RESULTS (CENTRAL BH) 

Note the typographical error in Metzger & Fernandez (2014) 
that reverses the order of the timescales. 


In this section, we present the results for the simula¬ 
tions where it is assumed that a black hole forms imme¬ 
diately upon merger and is present in every snapshot. 
Table 3 lists the times at which we simulate MC neu¬ 
trino transport and the corresponding global fluid and 
neutrino radiation properties. In what follows, we will 
probe the neutrino radiation field and its interaction with 
the fluid, and try to explain differences between the Eull 
MC simulations, the Simple MC simulations, and the 
leakage data of ME 14. We do most of our comparisons 
with snapshots from a time of 3 ms after the start of the 
dynamical simulation as differences between the meth¬ 
ods are most striking then, though the composition and 
amount of ejecta are determined by the long-term evolu¬ 
tion. 

Mentioning one caveat is in order. The fluid snap¬ 
shots were evolved in the dynamical simulations using the 
leakage neutrino treatment of ME 14, and we simply per¬ 
form our MC transport on snapshots of the evolved fluid. 
Though the Simple MC transport employs the same set 
of neutrino interactions, the geometry and spectral shape 
of the neutrino radiation field are different from what is 
assumed in the leakage scheme. Eluid may be near ther¬ 
mal and weak equilibrium with neutrinos in the leakage 
scheme, but this is not necessarily true after switching 
to MC transport. This potential discrepancy can easily 
cause artificially large or small heating and leptonization 
rates. If the fluid were evolved with an MC treatment 
instead, it would likely be much closer to equilibrium 
with the MC neutrinos, and the rates might not be as 
high. Because of this, our comparisons between leak¬ 
age and MC results indicate the qualitative effects, such 
as faster cooling and leptonization rates, but the magni¬ 
tudes of the differences are likely not reliable. The effects 
of MC transport on the end results of dynamical simu¬ 
lations are thus difficult to determine. The dynamical 
simulations also begin with a disk of uniform electron 
fraction Yg = 0-1, which is not initially in equilibrium 
with either leakage or MC neutrinos. Addressing this 
out-of-equilibrium issue requires the dynamical simula¬ 
tions to begin before the merger and to be coupled to 
MC neutrino transport, which we leave to future work. 

3.1. Neutrino Radiation Field (Central BH) 

In Eigure 2 we show the spatial distribution of the neu¬ 
trino radiation field at t = 3 ms for both Simple and 
Eull neutrino physics. Though the plot includes all neu¬ 
trino species, the radiation is dominated by electron anti¬ 
neutrinos. Most of the neutrino energy comes from the 
inner regions of the disk close to the central object, and 
the dense disk casts a shadow that reduces the neutrino 
luminosity and energy density at large radii near the 
equator. The inner boundary also blocks neutrinos from 
moving to the other side of the disk and creates a polar 
shadow. The elastic electron scattering in the Eull simu¬ 
lations results in higher opacities, which in turn deepens 
the equatorial and polar shadows. The Lorentz trans¬ 
formation of neutrinos in fluid moving at around 0.6c 
near the inner boundary increases the average energy of 
neutrinos emitted from the hot inner disk by ^ 30%. 
Additionally, the neutrinos are beamed in the azimuthal 
direction, causing fewer of the higher-energy neutrinos 
coming from the inner parts of the disk to be present 
in the polar regions and more to be present along the 




Figure 2. Neutrino energy density and average energy 
(Central BH) at t = 3 ms. Left hemisphere: Neutrino energy 
density, summed over all species and multiplied by to remove 
effects of distance from the center. Right hemisphere: energy 
density-weighted average neutrino energy, averaged over all species. 
Bottom hemisphere: Simple MC results. Top hemisphere: Full 
MC results. The black curve is the p = 10® g cm“^ contour, be¬ 
low which Sedonu opacities and emissivities are set to zero. The 
outer radius on the plot is at 250 km and the inner radius is 30 km. 
The neutrino radiation field is very asymmetric and sensitive to 
the included physics. The disk casts a shadow as higher-energy 
neutrinos are preferentially absorbed. Much more asymmetry is 
present when the Full suite of physics is included. 


45° radial. With either set of physics, this is different 
from the neutrino radiation field described by Fernandez 
& Metzger (2013), which becomes spherical at large dis¬ 
tances. 

Profiles of the neutrino radiation field split into the dif¬ 
ferent neutrino species are shown in Figure 3. There is 
some noise at small radii due to a relatively small num¬ 
ber of simulated neutrino packets present there. Along 
the equator, normalized neutrino energy density and av¬ 
erage energy are much higher close to the black hole than 
farther out in the disk, since the higher-energy neutrinos 
are preferentially absorbed by the disk. Moving radially 
along the pole, the energy density increases quickly as 
more of the disk becomes visible. However, the aver¬ 
age energy is always close to the average energy escaping 
from the disk listed in Table 3. In all directions, the neu¬ 
trino species follow the same hierarchy, such that electron 
anti-neutrinos everywhere contribute most to the energy 
density and heavy lepton neutrinos everywhere have the 
highest average energy. 

^emit: ^^escape: {E z^jemit): S^nd (-E'zy,escape) 1^ Table 3 dc 
scribe the global lab-frame properties of the neutrinos 
that are emitted and of those that escape through the 
outer boundary. In the snapshot at 3 ms, there is more 
energy emitted as electron neutrinos than as electron 
anti-neutrinos, but this is before the initial data in the 
simulations of MF14 has had any time to come into a 
quasi-equilibrium with the viscous heating and neutrino 
interactions. At all other times, electron anti-neutrino 
emission is stronger, reflecting a tendency of the fluid to 
relax to a higher electron fraction. The subsequent emis¬ 
sion rates of all species decrease with time as the disk 
loses mass and cools. Heavy lepton neutrinos interact 
much more weakly with the fluid than do electron neutri- 



Figure 3. Neutrino Radiation Profile (Central BH) from 
the Full MC simulation at t = 3 ms for all three simulated neutrino 
species. The neutrino radiation field is asymmetric and dominated 
by electron anti-neutrinos. Top: neutrino energy density along 
the pole (solid lines) and the equator (dashed lines), multiplied by 
to remove effects of distance from the center. Bottom: energy 
density-weighted average neutrino energy along radial lines. The 
green Ux curves represent the sum of all four heavy lepton neutrino 
species. 



Figure 4. Neutrino Spectra (Central BH) from the Full MC 
simulation at t = 3 ms. Dashed lines are spectra of each neutrino 
species escaping from within 10° of the equator, while the solid 
lines are those from within 10° of the 45° cones, normalized by 
the solid angle covered by the respective regions. Overplotted for 
both directions (distinguished by proximity to the data curves) are 
dotted zero-chemical potential blackbody curves with the same to¬ 
tal ffux and average energy as the measured spectrum. The large 
dot on the blackbody curve indicates this average energy. For 
smoothness, the spectra are taken from the 2xEnergy run in Ap¬ 
pendix A. The escaping neutrino radiation is somewhat nonthermal 
and asymmetric. 
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nos and anti-neutrinos both in absorption and emission 
since they participate only in neutral-current reactions. 
This, combined with the low optical depths that prevent 
a blackbody distribution from building up, causes the 
heavy lepton neutrinos to always be subordinate to elec¬ 
tron neutrinos and anti-neutrinos in energy density and 
in fluid heating and cooling. 

If we assume the neutrinos form a zero-chemical poten¬ 
tial blackbody distribution as is done in MF14, we can 
relate temperature to average energy through 


{Eu) 


J^E^BeMT) dE, 

f^BEA0,T)dE, 


A.llkbT , 


(13) 


where Ej^ is the neutrino energy, is the Boltzmann 
constant, and is the neutrino blackbody func¬ 

tion at temperature T and chemical potential ijl (Equa¬ 
tion Bl). While the density-weighted average fluid tem¬ 
perature at t = 3 ms is around 3 MeV, the average emit¬ 
ted neutrino energy is between 26 and 29 MeV for all 
species. Thus, most of the neutrinos are created in the 
hottest regions of the disk very close to the black hole. 
The opacity to neutrinos scales approximately like E‘1 
(e.g.. Burrows et al. 2006) causing more higher-energy 
neutrinos to be absorbed and the average energy of escap¬ 
ing neutrinos to be much smaller than that of the emitted 
neutrinos. The heavy lepton neutrinos have the coolest 
emission temperature but the hottest escape tempera¬ 
ture, since electron neutrinos and anti-neutrinos have 
much larger opacities and higher-energy neutrinos are 
preferentially absorbed. 

The disk’s self-irradiation in the leakage scheme is cal¬ 
culated as in MF14, and the global properties are sum¬ 
marized in Table 1. The temperature of the radiation 
in the leakage scheme is determined by an emissivity- 
weighted average and is the same for both electron neu¬ 
trinos and anti-neutrinos. The average energy in Table 1 
is computed from a zero-chemical potential blackbody 
of this temperature using Equation 13. Both the average 
energy and volume-integrated emission luminosities from 
the leakage data are much lower than those computed by 
Sedonu, since the “emission” in the leakage scheme ap¬ 
proximately accounts for immediate re-absorption in the 
same grid cell. The leakage and MC emission quantities 
then do not represent the same physics, but the differ¬ 
ence further illustrates that the higher-energy neutrinos 
are re-absorbed locally while the lower-energy neutrinos 
are able to escape. 

Using instead the Simple set of physics described in 
Table 2 does little to bring the leakage (Table 1) and 
Sedonu (Table 3) results closer together. However, it 
does result in significant deviations from the Full set of 
physics. In the following, we exclude individual pieces of 
physics to pinpoint the origin of the differences in the t = 
3 ms snapshot. The exclusion of scattering predictably 
does nothing to the properties of the created neutrinos in 
the disk, but by decreasing the optical depth, allows more 
of the higher-energy neutrinos to escape. Ignoring special 
relativity results in a decrease of the emission luminosity 
and average energy of all species. Neglecting to correct 
for weak magnetism and recoil effects causes the electron 
anti-neutrino emission rate to increase by about 20%, but 
since most of the opacity is also increased by a similar 
amount, the escaping luminosity increases by only about 


3%. The other species are minimally affected. 

The low electron fraction throughout the disk causes 
electron neutrinos to be very likely to absorb onto neu¬ 
trons, allowing few of them to escape. Their escape lumi¬ 
nosity in Table 3 is around an order of magnitude lower 
than their emitted luminosity, indicating an average op¬ 
tical depth of r ^ 3 for electron neutrinos, r ^ 1 for elec¬ 
tron anti-neutrinos, and r ^ 1 for heavy lepton species. 
Figure 4 shows neutrino spectra that further demonstrate 
the asymmetry of the escaping neutrino radiation. The 
leakage scheme does not yield data with which we can 
directly compare our escape spectra. However, we show 
that the spectra appear qualitatively similar in shape to 
zero-chemical potential blackbody spectra with the same 
average energy and total flux (dotted lines in Figure 4), 
though the MC spectra are somewhat pinched with peak 
energies higher by a few MeV. 

3.2. Neutrino-Fluid Interaction (Central BH) 

Different treatments of neutrino effects can have a sig¬ 
nificant impact on the fluid evolution. Panel A of Fig¬ 
ure 5 shows rapid cooling in the densest parts of the disk 
and net heating above and below the equatorial plane in 
both the MC and leakage results at t = 3 ms. However, 
MC transport results in faster heating directly above the 
disk by more than an order of magnitude, a smaller cool¬ 
ing region, and much faster cooling on the equator than 
leakage, making the leakage heating nearly invisible in 
Panel A of Figure 5. The differences are largely due to 
the approximate nature of the disk self-irradiation and 
leakage scheme of MF14, the accuracy of which suffers 
especially at the midplane of the disk. For efficiency, the 
leakage scheme calculates the optical depth at a given 
point to he r = K min{r, i^_L, ^||}, where hz is the opac¬ 
ity, r is the radius, and H± and are the vertical and 
horizontal pressure scale heights, respectively. This opti¬ 
cal depth calculation is only accurate to within a factor 
of a few. MC transport allows neutrinos to escape in 
any direction rather than just vertically and radially, in¬ 
creasing the amount of escaping neutrinos. Using Full 
neutrino physics dramatically increases the heating rate 
just above the hottest part of the disk due to special 
relativistic effects boosting the luminosity and average 
energy of neutrinos, and the larger global heating rate is 
also reflected in a smaller value of Cjy — Hu in Table 3. 

Viscosity, neutrino heating, and neutrino cooling all 
affect the thermal evolution of the disk. The relative im¬ 
portance of neutrinos and viscosity can be seen in the 
amount of mass for which neutrino heating is larger than 
viscous heating in Table 3. At the 3 ms Full MC snap¬ 
shot, 3.15 X 10is heated more strongly by neutrinos 
than viscosity, though after this time the number drops 
very quickly. Simple neutrino physics causes this mass 
to be 60% smaller, and in the Leakage simulation this 
mass is almost zero. In Panel B of Figure 5 we show 
Re = {l/e)de/dt^ the relative rate of change of internal 
energy including the viscous heating rate used in the dy¬ 
namical simulations of MF14, scaled by each point’s cur¬ 
rent energy density. The extra heated wings above the 
disk in the MC simulation cause the internal energy to 
change several times faster than leakage would suggest. 
More dramatically, the MC simulations show neutrino 
cooling dominating viscous heating along the equator, 
while the opposite is true in the leakage results. From 
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Figure 5. Neutrino-Fluid Interaction (Central BH) a,t t = Sms. In each plot, the right half shows results calculated with neutrino 
leakage in the dynamical simulations of MF14, while the left half is calculated by Sedonu. Each quadrant of Sedonu results depicts only half 
of the simulation domain. The top left quadrant uses the Full set of physics, while the bottom left uses the Simple set of physics. Panel A 
shows the difference between absorptive heating and emissive cooling. Panels B and C depict the relative rate of change of internal energy 
(including viscous heating) divided by internal energy and of electron fraction, respectively. Red represents a large positive rate of change 
while blue represents a large negative rate of change. Any rate of change whose magnitude is smaller than ls“^ is plotted as 0. The outer 
radius on each plot is at 250 km and the inner radius is 30 km. The black curve is the p = 10® g cm“^ contour, below which Sedonu opacities 
and emissivities are set to zero. Using MC neutrino transport would likely significantly affect the thermal and compositional evolution of 
the disk. 



Figure 6. Equilibrium electron fraction (Central BH) at 

t = 3 ms. Left Hemisphere: equilibrium electron fraction at which 
the net lepton number absorption rate is equal to the net lepton 
number emission rate. The equilibrium solver is unreliable below ~ 
1.5 MeV as the energy grid ceases to be able to resolve the neutrino 
distributions, so Te^goived “ at locations with a temperature less 
than this is plotted as zero. Right Hemisphere: rate of change 
of electron fraction, as in Figure 5, but separately depicting that 
caused by emission (top) and absorption (bottom). The black curve 
is the p = 10® g cm ^ contour, below which Sedonu opacities and 
emissivities are set to zero. The outer radius on the plot is at 
250 km and the inner radius is 30 km. 


this we would expect a stronger neutrino-driven wind, a 
thinner disk, and much faster disk cooling, though dy¬ 
namical simulations would be required to investigate this 
quantitatively. 

In a similar manner, we show Ry^ = dYejdt in Panel C 
of Figure 5. There is very little difference between the 
Simple and Full MC simulations, though both represent a 
significant departure from the leakage data. In all cases, 
electron fraction is increasing near the equator as the 
low-electron fraction fluid is emitting primarily electron 
anti-neutrinos. Above the equator, there is a pattern of 
regions of both increasing and decreasing electron frac¬ 
tion. The right half of Figure 6 shows that this pattern of 
increasing and decreasing electron fraction is caused by 
neutrino emission rather than absorption, indicating that 
variations in density and temperature cause the equilib¬ 
rium electron fraction to vary. In the left half of Figure 6, 


we show the difference between the current electron frac¬ 
tion and the equilibrium electron fraction (including neu¬ 
trino interactions, see Section 2). The initial conditions 
of the dynamical simulation began with an electron frac¬ 
tion of Ye = 0-1, leaving the center of the disk far from 
equilibrium, and the slower rates in the leakage scheme 
prevent it from coming into equilibrium more quickly. 

As the disk spreads, cools, and accretes onto the cen¬ 
tral BH, neutrinos affect the evolution of the disk and 
outflow more slowly. However, MC transport differs sig¬ 
nificantly from leakage for at least several tens of mil¬ 
liseconds. The volume-integrated neutrino cooling mi¬ 
nus heating through m the 3 ms snapshot is up to ^ 38% 
larger in the Full MC simulations than in the Leakage 
ones. At all later times, though, the Leakage simulations 
cool faster by up to 22% a similar percentage . Through¬ 
out the disk’s evolution MC results in a higher leptoniza- 
tion rate, up to 7.7 times that of the Leakage simulation 
at 3 ms. This is due to the way the leakage scheme treats 
optical depths, as discussed in Section 6. 

In Figure 7, we show the difference between the rates 
of change of internal energy due only to neutrinos {IZ'J 
and electron fraction (T^y^) calculated using leakage and 
using MC transport. Though the differences are most 
striking at 3 ms, we still see a significantly larger 7^' near 
the 45° radials and a larger near the midplane in the 
MC results at 30 ms. In these simulations, Sedonu takes 
the opacity at densities lower than 10^ g cm“^ to be zero, 
so the leakage scheme shows more neutrino heating for 
a very small amount of mass outside the disk. However, 
the left and right halves of Panel B of Figure 5 outside 
of the region covered by Sedonu appear almost identi¬ 
cal because in this region viscous heating is completely 
dominant. Comoving frame viscous heating is identical 
in all cases, but time dilation slightly modifies the heating 
rates in all but the Simple and NoRel cases. At 300 ms 
neutrino cooling by leakage is more efficient than Monte 
Carlo cooling, though the differences are only apparent 
at the densest part of the disk and are anyway dominated 
by viscous heating. 

For calculating opacities and interaction rates, NuLib 
requires as input an EOS to determine the chemical 
potentials of the particles involved in each interaction. 
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Figure 7. Leakage-MC Difference (Central BH). Full MC 
transport results differ significantly from leakage results for at least 
several tens of milliseconds. Difference between the Full and Leak¬ 
age relative rate of change of internal energy ignoring viscosity 
(top panel) and electron fraction (bottom panel), multiplied by 
the snapshot time in order to estimate the potential impact of im¬ 
proved neutrino transport on dynamical simulations. The black 
curve is the p = 10® g cm“^ contour, below which Sedonu opacities 
and emissivities are set to zero. The outer radius on the plot is at 
250 km and the inner radius is 30 km. In both plots, values larger 
than unity imply that if such rates continued for a time equal to 
the snapshot time, the difference between MC and leakage would 
be dynamically important. For this plot, data is taken from the 
2xEnergy run in Appendix A for increased solution accuracy. 



Thus, different EOS result in different interaction rates. 
In addition to the fiducial Helmholtz EOS (Timmes & 
Swesty 2000), we repeat the calculations using the LS220 
(Lattimer & Swesty 1991) and the H. Shen (Shen et ah 
2011) EOS. We find that the choice of EOS has no signifi¬ 
cant effect, as is demonstrated by the results summarized 
in Table 3. This is reassuring, since all of the neutrino 
emission and absorption occurs at sub-nuclear densities 
where the details of the treatment of the strong force are 
less significant. 

4. RESULTS (CENTRAL HMNS) 

Eollowing the merger of two neutron stars, the central 
object that forms may be a black hole, a stable neutron 



Figure 8. Neutrino energy density and average energy 
(Central HMNS) at t = 3 ms. The quantities shown are the same 
as in Figure 2, but using the HMNS snapshot. Left hemisphere: 
Neutrino energy density, summed over all species and multiplied by 
to remove effects of distance from the center. Right hemisphere: 
Neutrino energy density-weighted average energy, averaged over 
all species. Bottom hemisphere: Simple MC results. Top hemi¬ 
sphere: Full MC results. The black curve is the p = 10® g cm“^ 
contour, below which Sedonu opacities and emissivities are set to 
zero. The outer radius on the plot is at 250 km and the inner 
radius is 30 km. The neutrino radiation field is very asymmetric 
and sensitive to the included physics. The disk casts a shadow as 
higher-energy neutrinos are preferentially absorbed. Much more 
asymmetry is present when the Full suite of physics is included. 
Both the energy density and average neutrino energy are higher 
than in the BH case. At this point in time, the luminosities and 
average energies of from the central HMNS are set to 

{10.2,10.2, 40.8} B s-i and {16.4, 20.5, 20.5} MeV, respectively. 


star, or an only temporarily-stable hypermassive neutron 
star (HMNS), depending on details of the equation of 
state and the object’s mass (e.g., Kaplan et al. 2014). In 
this section, we bracket the parameter space by repeat¬ 
ing the analysis of the previous section, but with simu¬ 
lations including an HMNS assumed to be permanently 
stable. The inner boundary, which models an HMNS by 
reflecting matter, prevents mass from accreting through 
it and leads to a disk that stays hot and massive for a 
much longer time. As did ME14, we assume the HMNS 
emits neutrinos with a zero-chemical potential blackbody 
distribution and the average energies listed in Table 1. 
Heavy lepton neutrinos, when present, have the same lu¬ 
minosity and average energy as electron anti-neutrinos. 
Table 3 lists the simulations we run and the correspond¬ 
ing global properties of the fluid and radiation field. 

4.1. Neutrino Radiation Field (Central HMNS) 

As seen in Eigure 8, in the presence of an HMNS the 
neutrino radiation field at t = 3 ms shows the same disk 
and polar shadows and relativistic beaming as when a 
BH is present, though the neutrino energy densities and 
average energies are somewhat higher due to the hot¬ 
ter fluid and extra irradiation from the HMNS. The Eull 
physics neutrino radiation field shows higher energy den¬ 
sities by a factor of ^ 1.5 in the free streaming regions 
outside the disk than the Simple physics neutrino radia¬ 
tion field. This is due mostly to the production of copious 
amounts of heavy lepton neutrinos in the HMNS, which 
have comparatively small cross sections and so are able 
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to pass much more easily through the disk. As iu the BH 
case, special relativity iucreases the average ueutriuo eu- 
ergy by up to ^ 30%, especially 45 degrees from the pole. 

^emit: ^escape? )5 aud escape) 1^ Xable 3 de¬ 

scribe the global lab-frame properties of the ueutriuos 
iu the simulatious with au HMNS. Siuce the iuitial disk 
couditious are the same for both dyuamical simulatious, 
the properties of the ueutriuos emitted from the disk 
at t = 0 ms are also ideutical. The extra irradiatiou 
from the HMNS results iu a slower uet cooliug rate of 
the disk at early times. After the first 3 ms, the large 
disk mass arouud the HMNS causes the uet cooliug rate 
to be much larger, though the higher temperatures aud 
amouut of irradiatiou by uearly equal uumbers of electrou 
ueutriuos aud auti-ueutriuos imposed from the HMNS 
bouudary couditiou (see Sectiou 2.3.1) cause a slower 
chauge iu electrou fractiou thau iu the BH case. The 
high deusities aud viscously-amplified temperatures uear 
the HMNS cause 1.9 times as much ueutriuo euergy to be 
emitted from the disk, though the vast majority is imme¬ 
diately re-absorbed. A combiuatiou of the HMNS’s extra 
radiatiou aud the higher disk lumiuosity causes the es- 
capiug ueutriuo lumiuosity to be larger by factors of 2.8, 
1.5, aud 90 for z/g? aud respectively. The euer- 
gies of the escapiug ueutriuos are similar to those iu the 
BH case, but the Ux average euergy is decreased due to 
dilutiou from the HMNS. Note that iu simulatious that 
iuclude the HMNS, the HMNS is much hotter thau the 
disk (e.g., Dessart et al. 2009), but we parameterize the 
ueutriuos beiug emitted from the HMNS for cousisteucy 
with the dyuamical simulatious of MF14. 

Figure 9 shows a complicated iuteractiou betweeu radi¬ 
atiou from the HMNS aud that from the disk. Aloug the 
poles there is au iuitial dip iu iuteusity as ueutriuos are 
absorbed by a layer of matter just outside of the HMNS. 
Moviug outward aloug the pole, as the disk comes iuto 
view, the electrou ueutriuos aud auti-ueutriuos emitted 
from the disk bump their respective iuteusities up agaiu. 
Though they iuteract more weakly thau electrou ueutri¬ 
uos aud auti-ueutriuos, heavy leptou ueutriuos are scat¬ 
tered by the disk, as seeu iu the divergeuce of the pole 
aud equatorial euergy deusity. Iu fact, the scatteriug 
combiued with Doppler boostiug iu the Full simulatiou 
cause the average heavy leptou ueutriuo euergy peak iu 
the iuuer regious of the disk, aud the scattered ueutri¬ 
uos eveu cause radially iucreasiug average ueutriuo eu¬ 
ergy aloug the poles. Below 150 km aloug the equator, 
electrou aud heavy leptou ueutriuo iuteusities decliue as 
their emissiou from the HMNS is absorbed by the disk. 
The deuse part of the disk below 60 km is such a stroug 
emitter of electrou auti-ueutriuos, however, that the iu¬ 
teusity rises before falliug agaiu farther out. Iu all cases, 
the ueutriuo radiatiou field becomes free-streamiug after 
a radius of ^ 150 km, iudicated by horizoutal hues iu 
Figure 9. 

Figure 10 demoustrates the asymmetry of the ueutriuo 
radiatiou field aud its departure from a zero-chemical 
poteutial blackbody. At all times, the irradiatiou from 
the HMNS makes the uet escapiug flux from heavy lep¬ 
tou ueutriuos comparable to that from electrou auti- 
ueutriuos, iu coutrast to the BH case. 

Just as with the BH suapshots, usiug Full physics adds 
a scatteriug opacity which preveuts ueutriuos from es¬ 
capiug as easily as with Simple physics, causiug Sim- 



Radius [km] 

Figure 9. Neutrino Radiation Profile (Central HMNS) 

from the Full MC simulation at t = 3 ms for all three simulated 
neutrino species. The quantities shown are the same as in Fig¬ 
ure 3. The neutrino radiation field is asymmetric and dominated 
by electron anti-neutrinos. Top: neutrino energy density along 
the pole (solid lines) and the equator (dashed lines), multiplied by 
to remove effects of distance from the center. Bottom: energy 
density-weighted average neutrino energy along radial lines. The 
green Ux curves represent the sum of all four heavy lepton neutrino 
species. Note the difference in the y-axis scale compared with Fig¬ 
ure 3. The hierarchy between neutrino species is shuffled at small 
radii due to competing emission from the HMNS and the disk, and 
from disk absorption. 


pie physics to allow for a larger cooliug aud leptouiza- 
tiou rate, as well as uaturally larger escape lumiuosi- 
ties. Full physics also iucreases the ueutriuo creatiou 
rate through weak maguetism correctious aud Loreutz 
trausformatious. Eveu though a large uumber of heavy 
leptou ueutriuos are produced by the HMNS, Figure 9 
shows that they coutribute much less to disk heatiug thau 
electrou auti-ueutriuos do. This is because they deposit 
euergy iuto the fluid ouly through NuLib’s approximate 
treatmeut of iuverse Bremsstrahluug aud ueutriuo pair 
auuihilatiou. Excludiug heavy leptou ueutriuos ouly re¬ 
sults iu the disk cooliug 3% more quickly. Otherwise, 
the exclusiou of various elemeuts of physics has the same 
effect as iu the BH suapshots. 

4.2. Neutrino-Fluid Interaction (Central HMNS) 

Eigure 11 describes the iuteractiou of ueutriuos with 
the backgrouud fluid. Near the equator, the structures 
of the heatiug rates, i?e, aud are very similar to the 
BH case. Pauel A shows the local ueutriuo heatiug rates, 
the volume iutegrals of which are displayed iu Table 3. 
Eull MC, as iu the BH case, shows slower iutegrated ueu¬ 
triuo cooliug thau the Simple MC simulatiou (factor of 
^ 0.91) but much faster cooliug thau the Leakage simu¬ 
latiou (factor of ^ 2.2). The relative effects of ueutriuos 
aud viscosity cau be seeu iu the amouut of mass for which 
ueutriuo heatiug is larger thau viscous heatiug. This 
mass is ^ 26% larger iu the Eull MC simulatious thau iu 
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Figure 10. Neutrino Spectra (Central HMNS) from the 
Full MC simulation at t = 3 ms. The quantities shown are the 
same as in Figure 4, though note that the vertical axis differs. The 
escaping neutrino radiation is somewhat nonthermal and asym¬ 
metric. Dashed lines are spectra of each neutrino species escaping 
from within 10° of the equator, while the solid lines are those from 
within 10° of the 45° cones, normalized by the solid angle covered 
by the respective regions. Overplotted for both directions (distin¬ 
guished by proximity to the data curves) are dotted zero-chemical 
potential blackbody curves with the same total ffux and average 
energy as the measured spectrum. The large dot on the blackbody 
curve indicates this average energy. For smoothness, the spectra 
are taken from the 2xEnergy run in Appendix A. The heavy lepton 
neutrinos are near blackbodies since most come from the central 
HMNS. 

the Simple MC simulation and ^5.6 times larger than 
in the Leakage simulation. This is visible in the relative 
sizes of the neutrino-heated regions above and below the 
disk. 

The rates of change of internal energy in Panel B of 
Figure 11 demonstrate that both Full and Simple MC 
neutrinos escape from and pass through the densest re¬ 
gions of the disk more easily than leakage allows, causing 
visibly faster heating near the equator beyond the disk. 
Similar to the heating rates, the volume-integrated lep- 
tonization rate indicated by the Full MC simulation is 
faster than that predicted by the Leakage simulation. 
The slight difference between the Simple and Full MC 
simulations is also reflected in Table 3. Simple physics 
decreases opacities and results in a faster increase in elec¬ 
tron fraction, mostly near the HMNS above the disk. 
Monte Carlo allows neutrinos to escape easier than they 
could in the Leakage simulation and hence has higher 
cooling and leptonization rates. 

Although the differences stem from the same effects, 
the Full MC net neutrino cooling rate is at times a fac¬ 
tor of 8 higher than the leakage net neutrino cooling rate 
(at t = 30 ms), though this is likely artificially high due to 
the out-of-equilibrium effects discussed in Section 3. Un¬ 
like in the BH case, MC net neutrino cooling minus heat¬ 
ing is larger than that from leakage calculations through 
300 ms since the disk retains its mass, and the greater 
ease of escape for MC neutrinos allows the disk to ab¬ 
sorb less energy from the HMNS neutrinos and to cool 
more quickly. In addition, the HMNS emits nearly equal 
numbers of electron neutrinos and anti-neutrinos, but the 
easier escape allowed to MC neutrinos means the HMNS 
neutrinos are not as effective at bringing the electron 


fraction up. At t = 3 ms, MC leads to faster cooling 
than leakage does by a factor of ^ 2.2 and leptonization 
by a factor of ^ 2.7. However, unlike the BH case, after 
the 3 ms snapshot MC leptonization is slower than leak¬ 
age. By the 3 s snapshot, the dynamical simulation has 
overshot MC’s equilibrium and so the MC disk is slowly 
heating rather than cooling. 

The disk mass remains in the HMNS simulation for al¬ 
most a hundred times as long as it does in the BH snap¬ 
shots, and neutrinos from both disk and HMNS emission 
play an important role for at least ten times as long as 
in the BH snapshots. Figure 12 compares the leptoniza¬ 
tion rates and the difference between cooling and heating 
rates due only to neutrinos through the 3 s snapshot. The 
data in the 3 ms quadrants effectively replicates informa¬ 
tion conveyed in Figure 11 by showing that MC transport 
allows the disk to cool faster, heats the regions above the 
disk faster, and allows the disk Tg lo change more quickly. 
In the 30 and 300 ms quadrants, leakage predicts that 
neutrinos are unable to cool the matter in the equatorial 
plane, but are able to cool the disk above and below the 
equator. MC transport, on the other hand, predicts some 
cooling in isolated domains on the equator and next to 
the HMNS at high latitudes, but neutrino heating bal¬ 
ances neutrino cooling in the mid-latitude disk regions. 
In addition, leakage predicts much stronger heating of 
the low-density polar regions, a trend which continues 
to be evident at later times. In the 3 s quadrants, very 
little mass remains in the disk and the predicted leakage 
rates are far in excess of the MC ones, just as they were 
in the low-density regions in previous snapshots. This is 
all consistent with the above statement that MC allows 
neutrinos to escape more easily than leakage does. 

We solve for the equilibrium electron fraction in Fig¬ 
ure 13. The results are very similar to the BH case. 
The leakage neutrinos effectively interact more strongly 
than MC neutrinos do, which causes the fluid in the low- 
density polar regions to an increase in electron fraction 
through absorption of electron neutrinos from the HMNS 
in the dynamical simulation, bringing it closer to equilib¬ 
rium. In the main region of the disk, the fluid is below 
the equilibrium electron fraction since there is a suffi¬ 
ciently large amount of mass that neutrinos have not yet 
been able to significantly raise the electron fraction. If 
significant neutrino processing occurs before the disk is 
formed (i.e., before the initial conditions of the dynami¬ 
cal simulations of MF14), the electron fraction would be 
higher and not as far from equilibrium. 

5. NEUTRINO PAIR ANNIHILATION 

We calculate neutrino pair annihilation rates in a post¬ 
processing step after neutrinos have finished propagating 
through the disk. The resulting rates are plotted in Fig¬ 
ure 14 for the 3 ms snapshots with both BH and HMNS 
backgrounds, and volume-integrated rates are given for 
every snapshot in Table 3. The NoPair simulations do 
not include pair processes in the NuLib tables, but the 
annihilation post-processing requires only neutrino dis¬ 
tribution functions and does not rely on the NuLib ta¬ 
bles. In the sparse polar regions, the density is low 
enough that annihilation would rapidly increase the tem¬ 
perature and entropy, which has the potential to gener¬ 
ate a rapid outflow. However, annihilation accounts for 
at most ^ 4% of the global energy gain/loss rate in any 
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Figure 11. Neutrino-Fluid Interaction (Central HMNS) at t = Sms. The quantities shown are the same as in Figure 5. In each 
plot, the right half shows results calculated with neutrino leakage in the dynamical simulations of MF14, while the left half is calculated 
by Sedonu. Each quadrant of Sedonu results depicts only half of the simulation domain. The top left quadrant uses the Full set of physics, 
while the bottom left uses the Simple set of physics. Panel A shows the difference between absorptive heating and emissive cooling. Panels 
B and C depict the relative rate of change of internal energy (including viscous heating) divided by internal energy and of electron fraction, 
respectively. Red represents a large positive rate of change while blue represents a large negative rate of change. Any rate of change whose 
magnitude is smaller than ls“^ is plotted as 0. The outer radius on each plot is at 250 km and the inner radius is 30 km. The black curve 
is the p = 10® g cm“^ contour, below which Sedonu opacities and emissivities are set to zero. Using MC neutrino transport would likely 
significantly affect the thermal and compositional evolution of the disk. 


snapshot with either central object, and so will not dra¬ 
matically affect the dynamics of most of the disk mass. 
In order to better estimate the role annihilation would 
have in driving a relativistic jet, we integrate only over 
regions less than 45° from the poles. This excludes an¬ 
nihilation in the bulk of the disk and in regions far from 
the poles, which are either too dense or too far from the 
pole to contribute to acceleration along the poles. The 
ratio of the total annihilation to that just within 45° of 
the poles can be as high as 120 in the BH snapshots 
(at 3ms) and 220 in the HMNS snapshots (at 30ms). 
Determining how much mass and energy is driven by an¬ 
nihilation, and whether this can actually produce a jet, 
would require including annihilation rates in the dynam¬ 
ical simulation. 

An order-of-magnitude estimate of the total energy de¬ 
posited in polar regions from neutrino pair annihilation 
can be found by time interpolating the volume-integrated 
values of annihilation rate given in Table 3 and integrat¬ 
ing assuming that 1-Lyy = where k and C are param¬ 
eters set to create a piecewise-continuous interpolation. 
Integrating this interpolation for the Full physics simula¬ 
tions, we find that the total amount of energy deposited is 
Euu,net = 2.2 X 10^^ erg after 300 ms for the BH case and 
Euu,net = 1-8 X 10^^ erg after 3 s for the HMNS case. If 
for the reasons above we include only the volume within 
45° of the poles as in Table 3, the integrated deposited 
energy becomes = 2.8 x 10^^ erg for the BH case 

and = 1-9 x 10^^ erg for the HMNS case. How¬ 

ever, 10^^ — 10^^ ergs“^ is required to launch a GRB jet 
(e.g., Lee & Ramirez-Ruiz 2007). Given the assumptions 
used in both the MG neutrino transport and the dynam¬ 
ical simulations of MF14, this calculation is certainly not 
accurate enough to definitively rule out the possibility of 
an annihilation-driven e+ — e~ jet, but it is on the lower 
end of this energy requirement. 

The effects of various approximations on the annihi¬ 
lation rate at 3 ms for both the BH and HMNS cases 
are also summarized in Table 3. Though special relativ¬ 
ity increases the average neutrino energy, it also beams 
the neutrinos along similar trajectories in the azimuthal 
direction, which decreases the relative angle between 


neutrinos and causes the NoRel annihilation rate to be 
higher than the Full rate by ^ 44% in the BH snap¬ 
shot and ^ 13% in the HMNS snapshot. Scattering off 
of rapidly moving fiuid boosts neutrino energies and so 
the annihilation rate in the NoScat simulation of the BH 
snapshot is ^ 5% lower than in the Full simulation. Ad¬ 
ditionally, scattering causes neutrinos that would other¬ 
wise have passed straight through the disk to be defiected 
up toward the polar regions, and together with the in¬ 
creased neutrino energy this results in a ^ 18% increase 
in the annihilation rate in the HMNS snapshot. The 
weak magnetism correction results in a smaller escape lu¬ 
minosity of electron anti-neutrinos, which decreases the 
annihilation rate by ^ 5% in the BH snapshot and ^ 2% 
in the HMNS snapshot. In the BH snapshots, so few 
heavy lepton neutrinos are produced that they provide 
essentially no additional contribution to pair annihila¬ 
tion rates, but the heavy lepton neutrinos emitted from 
an HMNS can increase the global annihilation rate by 
^ 20%. These numbers apply only to the 3 ms snap¬ 
shots, but indicate the direction and approximate rela¬ 
tive magnitude of the effect each piece of physics has on 
the instantaneous annihilation rate. 

6. DISGUSSION 

Alhough MG neutrino transport is a less approximate 
treatment of neutrinos than leakage, it should be noted 
that there are still many approximations being made. 
The largest is the neglect of general relativity, which 
would red/blueshift neutrinos moving outward/inward, 
respectively, and would bend the neutrino trajectories 
along geodesics. Previous works have indicated that this 
following of geodesics causes neutrino trajectories to in¬ 
tersect at higher angles, increasing the annihilation rates 
by at most a factor of two (Asano & Fukuyama 2000, 
2001; Miller et al. 2003; Birkl et al. 2007; Harikae et al. 
2010). Moreover, the step between depositing neutrino 
energy and building a jet cannot be determined by our 
stationary simulations. Neutrino annihilation in the very 
sparse regions would also result in a very large entropy 
per baryon and thus a very efficient r-process even in 
matter that is barely neutron-rich, but the amount of 
mass in the polar regions is so small that the amount 
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Figure 12. Leakage-MC Difference (Central HMNS). The 

quantities shown are the same as in Figure 7. Difference between 
the Full and Leakage relative rate of change of internal energy 
ignoring viscosity (top panel) and electron fraction (bottom panel), 
multiplied by the snapshot time in order to estimate the potential 
impact of improved neutrino transport on dynamical simulations. 
The black curve is the p = 10® g cm“® contour, below which Sedonn 
opacities and emissivities are set to zero. At 300 ms everything in 
the image is above the density cutoff. The outer radius on the plot 
is at 250 km and the inner radius is 30 km. In both plots, values 
larger than unity imply that if such rates continued for a time 
equal to the snapshot time, the difference between MC and leakage 
would be dynamically important. Full MC transport results differ 
significantly from leakage results at all times. For this plot, data is 
taken from the 2xEnergy run in Appendix A for increased solution 
accuracy. 


of r-process elements would be insignificant (Fernandez 
et al. 2015a). 

Second, the annihilation kernels we use do not account 
for final-state electron and positron blocking, and are 
valid only for neutrinos with energies much larger than 
the electron rest mass (see Appendix C). Third, neutri¬ 
nos are fermions, and Pauli exclusion in regions where 
neutrinos are degenerate should affect their trajectories 
(e.g., Janka et al. 1992). We account for the fermionic 
nature of neutrinos only in their interactions with mat¬ 
ter and not in their propagation, but the very low de¬ 
generacy of the neutrinos makes this a good approxima¬ 
tion. Fourth, scattering kernels are actually inelastic and 



Figure 13. Equilibrium Electron Fraction (Central 
HMNS) at f = 3 ms. The dense region of the disk is far from 
equilibrium, as is the ffuid in contact with the HMNS at high lat¬ 
itudes. The quantities shown are the same as in Figure 6. Left 
Hemisphere: equilibrium electron fraction at which the net lepton 
number absorption rate is equal to the net lepton number emis¬ 
sion rate. The equilibrium solver is unreliable below 1 MeV as the 
energy grid ceases to be able to resolve the neutrino distributions, 
so Tg^soived “ at locations with a temperature less than this 
is plotted as zero. Right Hemisphere: rate of change of electron 
fraction, as in Figure 5, but separately depicting that caused by 
emission (top) and absorption (bottom). The black curve is the 
p = 10® g cm“^ contour, below which Sedonn opacities and emis¬ 
sivities are set to zero. The outer radius on the plot is at 250 km 
and the inner radius is 30 km. 


anisotropic, but we treat them as elastic and isotropic, 
and include a correction factor to approximate “effec¬ 
tive” anisotropic scattering (e.g.. Burrows et al. 2006), 
and we ignore inelasticity for the sake of simplicity. Since 
scattering opacities can have a significant impact on the 
neutrino radiation field, proper treatment of inelastic 
and anisotropic scattering could become an important 
source of energy deposition as, e.g., in the context of 
core-collapse supernovae (Lentz et al. 2012). Finally, the 
opacities in the outskirts of the disk where T <^0.5 MeV 
depend on the composition, which is likely not in NSE. 
However, addressing these approximations is beyond the 
current capabilities of Sedonn and we defer to future 
work to evaluate their importance more carefully. 

We see very significant differences in the cooling and 
leptonization rates between MC transport and leakage. 
The rate of energy loss predicted by MC transport is con¬ 
sistently larger than that predicted by leakage, but some 
systematic error is not unexpected from such an approxi¬ 
mation and we expect the MC neutrinos to be out of equi¬ 
librium with the fluid evolved with the leakage scheme. 
In the BH snapshots, the global leptonization rates cal¬ 
culated by MC transport are ^ 7 times larger than those 
calculated by leakage through the 30 ms snapshot. Since 
neutrino opacities scale roughly with the square of the 
neutrino energy (e.g.. Burrows et al. 2006), the leakage 
scheme used by MF14 attempts to account for the energy 
loss rate by calculating the optical depth based on the 
opacity of the fluid at the neutrino energy 
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Figure 14. Neutrino Annihilation Rates at t = 3 ms for the 

BH case (top panel) and the HMNS case (bottom panel). Left 
hemisphere: annihilation rate per unit volume. Right hemisphere: 
the same annihilation rate per unit mass. Top hemisphere: Full 
MC simulations. Bottom hemisphere: Simple MC simulations. 
The outer curve is the p = 10® g cm“^ contour, below which Sedonu 
opacities and emissivities are set to zero. The outer radius on the 
plot is at 250 km and the inner radius is 30 km. Most annihilation 
occurs in the dense disk, but is most significant per unit mass along 
the poles. 


where J-nift) are Fermi integrals of order n Equation B3). 
Note that here we assume the chemical potential is zero, 
as do MF14. When applied to energy escape, this ac¬ 
counts for the fact that low-energy neutrinos are able to 
escape more easily than higher-energy neutrinos due to 
the scaling of the opacity with neutrino energy. However, 
this choice of energy is designed to properly account for 
energy loss, not lepton number change. If we repeat the 
same exercise, but replace By{T) with By{T)/Ey to rep¬ 
resent number escape rather than energy escape, we get 


{E.) = 


lfTE^B^{0,T)dE,/E, 


fTB,{0,T)dEJE, 

= yj'4(0)/J'2(0) = 3.59 keT . 


(15) 


Thus, the average energy to use when calculating opac¬ 
ities to account for lepton number loss is about 11% 


smaller than that used when accounting for energy loss 
rate. Since the opacity for lower energy neutrinos is 
lower, using the same mean opacity for number and en¬ 
ergy escape causes the leakage scheme to underestimate 
the number of escaping leptons. The leakage scheme 
could be made more consistent by calculating separate 
optical depths for neutrino energy and number escape 
(Ruffert et al. 1996) or by calculating separate optical 
depths for each energy bin (e.g., Perego et al. 2014 in 
the context of the isotropic diffusion source approxima¬ 
tion). Additionally, since this issue is a spectral effect 
rather than a geometric one, it could be accounted for in 
many of the more sophisticated energy-dependent trans¬ 
port schemes, such as spectral two-moment transport 
(e.g., Shibata et al. 2011; Cardall et al. 2013; Just et al. 
2015b). 

One of the main differences between the effects of neu¬ 
trinos calculated by leakage versus those calculated by 
MC transport is the amount of heating above the disk 
where densities are relatively low. ME 14 argued that 
neutrinos are unable to drive a significant wind, but can 
affect the composition, especially above the disk, which 
increases the electron fraction of the viscously driven 
ejecta. To estimate the potential impact of the increased 
heating in our MC simulations, we look at the amount of 
mass Ml, for which neutrino heating dominates over vis¬ 
cous heating. Mi, is listed in Tables 1 and 3. In the BH 
snapshots, the leakage results indicate that essentially no 
mass is heated more strongly by neutrinos than by vis¬ 
cosity. However, for the first several tens of milliseconds, 
MC transport shows neutrinos being dynamically impor¬ 
tant in 11-50% of the mass in same snapshots (though 
quickly approaching zero after that). In the HMNS snap¬ 
shots, Ml, increases with time according to leakage (^ 8% 
of the disk mass at 0 ms to ^ 73% at 3 s), but decreases 
according to MC transport (^ 42% at 3 ms to essentially 
none at 3 s). This is largely due to the disk spreading 
out in the HMNS simulations such that much of the disk 
mass is below the minimum density for which neutrino 
interactions are accounted for in Sedonu. Though essen¬ 
tially all of the mass is still above the minimum density 
at 30 ms, by 300 ms 75% of the disk mass is above the 
minimum density (average density is 10^^ g cm“^), and 
by 3 s only 0.09% is above the minimum density (average 
density is 10^g cm“^). 

There is still much to be done before predictive simu¬ 
lated kilonova light curves become available, but the dif¬ 
ferences we see between the leakage and MC results could 
have dramatic implications for the elements formed in 
the ejecta and the resulting light curve. Previous studies 
indicate that the production of heavy r-process elements 
requires electron fractions below Ye ^ 0.2 — 0.3 (e.g., 
Wanajo et al. 2014; Kasen et al. 2015). Our Figures 7 
and 12 suggest that significant increases to the electron 
fraction of the disk ejecta are possible with Monte Carlo 
neutrino transport, since it results in the matter outside 
of the disk being more strongly neutrino processed. A 
weak r-process would still make elements up to A ^ 90 
in electron fractions up to Ye ^ 0.4 if the entropy is suf¬ 
ficiently high (Wanajo et al. 2014). However, the lack 
of a strong r-process in the disk wind would imply a 
stronger early blue peak in the kilonova light curve if the 
merger is observed from a polar direction where the disk 
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is not obscured by the lanthanide-rich dynamical ejecta 
(MF14;Kasen et ah 2015). 

We can apply the interpolation and integration scheme 
we used in Section 5 to the annihilation rates listed in 
Table 1 of Dessart et al. (2009) (assuming no BH spin) 
to estimate the total annihilation energy deposited for 
the 100 ms simulation to be 5.1 x 10^^ erg. This is about 
two orders of magnitude smaller than our estimate in 
Section 5 using the entire domain in the HMNS case, 
but is very similar to the estimate using only the regions 
within 45° of the poles. However, a direct comparison is 
somewhat difficult for the following reasons. (1) Dessart 
et al. (2009) do not have an inner boundary condition, (2) 
their HMNS is ^ 0.5Mq less massive than the one we as¬ 
sume, (3) the HMNS luminosity is an order of magnitude 
more luminous at 30 ms, (4) the disk is ~ 7 times more 
massive (with correspondingly higher densities), (5) they 
neglect viscous heating, and (6) they use a density cut¬ 
off of p = lO^^gcm”^ for calculating annihilation rates 
rather than an angle from the poles. This cutoff serves 
to exclude the HMNS and dense inner disk from the an¬ 
nihilation calculations, since energy deposited there is 
effectively trapped and unable to contribute to outflows. 
Since our disk is so much less massive, all of our disk 
mass has density p < 10^^ gcm“^. Given the vast differ¬ 
ences in the background fluid with which the annihilation 
calculations were performed, the differences in the anni¬ 
hilation rates are reasonable and to be expected. 

The annihilation rate, however, depends on the prod¬ 
uct of the neutrino and anti-neutrino intensities, and 
so is sensitive to changes in the neutrino luminosity. 
Though Dessart et al. (2009) see neutrino luminosities 
similar to ours, other studies show somewhat higher lu¬ 
minosities of all species. Foucart et al. (2015) simulate 
a NS-BH merger and find electron neutrino luminosi¬ 
ties of < 100 B s~^, electron anti-neutrino luminosities 
of < 300 B s“^ and collective heavy-lepton neutrino lu¬ 
minosities of < 100 B s~^ for a few tens of milliseconds. 
Sekiguchi et al. (2015) simulate a NS merger including 
the HMNS and see see similar luminosities over a similar 
time. In both cases, the luminosities are larger than ours 
by a factor of a few, which has the potential to increase 
the annihilation rate by an order of magnitude. If these 
luminosities are closer to those in nature than those com¬ 
puted by us and Dessart et al. (2009) and the geometry 
of the emission favors increased annihilation rates, there 
may yet be hope for the neutrino annihilation-powered 
GRB model. 

Dynamical simulations with MG neutrino transport (or 
other methods more sophisticated than leakage) are re¬ 
quired to determine the true long-term effects of the in¬ 
creased cooling and leptonization rates. For instance, al¬ 
though our results show that MG results in faster global 
cooling at all times, what may happen in a full simula¬ 
tion is faster cooling at early times and slower cooling 
at late times, since the disk will have become cold much 
faster. However, MG transport is currently too compu¬ 
tationally expensive to be used at every timestep in a 
three-dimensional dynamical calculation. Other trans¬ 
port methods like energy-dependent two-moment trans¬ 
port (e.g., Thorne 1981; Shibata et al. 2011; Just et al. 
2015b) will be able to account for the spectral effects 
and many of the geometric ones, but in this approxi¬ 


mation one must choose an otherwise undetermined clo¬ 
sure relation to close the system of equations. Several 
physically-motivated analytic closures and variable Ed¬ 
dington factor methods (e.g., Gardall et al. 2013) have 
been proposed, but any method with a local closure (i.e., 
one that is determined only by the radiation in the cur¬ 
rent grid cell) introduces a nonlinearity into the transport 
equation that leads to unphysical radiation shocks (e.g., 
Olson et al. 2000). In the future, it may be possible to 
find a closure treatment that maintains the efficiency of 
the two-moment transport scheme while accurately re¬ 
producing Monte Garlo results. 

7. CONCLUSIONS 

We present the new open-source, steady-state, special 
relativistic Monte Carlo neutrino transport code Sedonu. 
It efficiently calculates the energy- and angle-dependent 
neutrino distribution function on multi-dimensional fluid 
backgrounds and solves for the equilibrium fluid temper¬ 
ature and electron fraction. We have simulated neutrino 
transport through snapshots of post-merger disks using 
Monte Carlo (MG) techniques with various elements of 
physics. We compare the results to the leakage scheme 
used in the original dynamical simulations by MF14. 
Since the Monte Carlo neutrinos are out of equilibrium 
with the fluid evolved with the leakage scheme in the 
dynamical calculations, we believe that the qualitative 
trends we indicate are robust. However, determining the 
magnitudes of the differences between the two methods 
would require MG transport to be coupled to the fluid 
evolution. In light of this, we summarize our findings 
below. 

1. Compared with leakage, MG transport results in 
global cooling and leptonization rates that are 
higher than those predicted by leakage during the 
optically-thick disk stage. If the disk is optically 
thin with a central BH, MG cooling is slower and 
leptonization is faster than leakage. If the disk is 
optically thin with a central HMNS, MG cooling 
is faster and leptonization is slower than leakage. 
This suggests a stronger blue component of the 
kilonova light curve if viewed from polar angles. 

2. MG exhibits up to an order of magnitude stronger 
neutrino heating above the disk that could increase 
the strength of a neutrino-driven wind. 

3. In the disk midplane, cooling via MG neutrinos 
dominates viscous heating through the 30 ms snap¬ 
shot with either central object, in constrast to the 
leakage results. 

4. The neutrino radiation field at large radii is very 
asymmetric, with most of the radiation escaping 
around 45° from the equator. 

5. The peak energies of the neutrino distribution func¬ 
tions can be shifted by a few MeV higher (z^e, ^e) oi* 
lower {vx) from the peak of a zero-chemical poten¬ 
tial blackbody with the same average energy and 
flux. 

6. Neutrino pair annihilation deposits an order of 
magnitude more energy with a central HMNS (^ 
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1.9 X 10^^ erg) than with a central BH 2.8 x 
10^6 erg), though this is still unlikely to be sufficient 
to drive a jet, though higher neutrino luminosities 
could make it plausible. 

7. Special relativity increases the average energies of 
escaping neutrinos by around 1 — 3 MeV and beams 
higher-energy neutrinos away from the poles. The 
inclusion of heavy lepton neutrinos, pair processes, 
scattering, weak magnetism, and variations in the 
equation of state have together at most a 10 % effect 
on the integrated cooling and leptonization rates. 
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APPENDIX 

A. MICROPHYSICS RESOLUTION STUDY 

To ensure that our results are robust against our choices of numerical discretization, we repeat transport simulations 
for the Sms snapshots with different discretizations, all using the Eull set of physics, and compare to the original 
simulations in Table 4. The discretization of the NuLib tables was tested by in turn doubling the neutrino energy, 
matter density, matter temperature, and electron fraction grid resolution. We also in turn double the angular resolution 
of the neutrino distribution functions in each cell, which only has an effect on the annihilation rates. Einally, we double 
the number of steps each neutrino packet must take by changing dceii to be 0.2 rather than 0.4 times the cell’s smallest 
dimension (see Section 2 for details). 

The differences between simulations with enhanced microphysics resolution in Table 4 and the originals in Table 3 
are all small, indicating that our discretization is sufficiently fine. Increasing the angular resolution of the annihilation 
kernels causes annihilation rates to drop slightly, supporting the supposition that most of the annihilation is from 
small incident angles. 

Increasing the number of Monte Carlo neutrino packets does not introduce any systematic change, but rather only 
reduces the size of random fluctuations in results when the same simulation is run multiple times. We use 2 — 4 x 10^ 
packets in order to keep random fluctuations of the results in Table 3 at ^ 0.1%. 

B. CODE TESTS 

We perform a pair of related tests to ensure that our transport and equilibrium finding methods arrive at the correct 
answer. In the first test, we demonstrate that we are able to reach a blackbody distribution function at high optical 
depth. In the second, we demonstrate that our equilibrium solver settles to the correct values of temperature and 
electron fraction. 


B.l. Blackbody Generation 

In this test we use a single unit-volume fluid cell with periodic boundary conditions. We give the cell a density, 
temperature, and electron fraction, and observe the resulting neutrino radiation that builds up within the cell, similar 
to Tubbs (1978). The neutrinos should settle into a Eermi blackbody distribution given by 




El/h^c^ 


(Bl) 


in CCS units of ergs“^cm“^MeV~^sr“^, where Ey is the neutrino energy (different from the packet energy Ep men¬ 
tioned in the main text), fi is the neutrino chemical potential, and T is the fluid temperature. This is identical to the 
photon blackbody, with two exceptions. The sign in the denominator originating from the Eermi-Dirac distribution, 
and there is no factor of 2 in the numerator because only left-handed neutrinos have been observed. Integrating over 
neutrino energy and angle gives the blackbody energy density 
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where 
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are Fermi integrals of order n. In the special case of /u = 0, this has a simple analytical solution akin to the Stefan- 
Boltzmann law, 
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Figure 15. Blackbody Generation Test: The equilibrium total neutrino energy density Ciy as a function of fluid density p, temperature 
T, and electron fraction Ye- The fluid variables are varied around p = 10^® g cm“^, T = 3.1623 MeV, and Ye = 0.3. The Sedonu-calculated 
energy densities (black crosses) match the theoretical ones (red lines) in a wide range of fluid conditions. Regions of mismatch at low 
and high temperatures occur as the peak of the neutrino distribution approaches the low and high neutrino energy limits of the NuLib tables. 


We let the fluid emit neutrinos and allow them to propagate until they are absorbed and compare the resulting 
neutrino energy density to Equation B2. Because scattering opacities become much larger than absorption opacities at 
high densities, we use the NoScat physics for efficiency. We expect the net energy density to be the sum of contributions 
from each of the 6 neutrino species, constrained by and = 0, where Vx represents any of the four 

heavy lepton neutrino/anti-neutrino species. The equilibrium can be taken directly from the EOS for a given 
{p, T, Ye}- Eigure 15 demonstrates the match between the predicted energy density and that calculated by Sedonu for 
many values in each direction of {p, T, Yg}- The plots appear very similar for all equations of state, but we use the 
Helmholtz EOS to complement the results plotted in the main text. The computed and theoretical results disagree 
at low and high temperatures, where the neutrino distribution functions extend past the energy limits in the NuLib 
tables. 


B. 2 . Blackbody Irradiation 

Rather than determining what radiation field is in equilibrium with the input fluid properties as in the previous test, 
we determine what fluid properties are in equilibrium with the input radiation. That is, we solve for the equilibrium 
properties of fluid that is allowed to relax in a bath of blackbody neutrinos. We set up a thin shell of fluid {dr/r = 10“^), 
apply a reflective outer boundary condition, and emit blackbody neutrinos from an absorbing inner boundary specified 
by a neutrino temperature Tinput and electron neutrino chemical potential input- The chemical potentials of all 
other neutrino species satisfy the constraints detailed in the previous test, namely that yUzyg,input = input and 

input = 0- Unlike in the main text, the luminosity of each species s is determined by the input temperature and 
chemical potential of the blackbody neutrinos being emitted from the inner boundary, such that 

Cs = dTT^r^ ,input 5 Tinput 5 

Eu,i 

where r is the radius of the inner boundary, is the center of energy bin i, and is the width of energy bin i. 

We then iteratively relax the fluid to its equilibrium temperature and electron fraction as described in the following. 
Each iteration is done as described in Section 2. If we are solving for temperature, we then set the temperature of 
each grid cell to + dAT^, where and are the temperatures for the iterations i + 1 and i, respectively, 

d = 0.3 is a somewhat arbitrary damping factor, and ATi is the difference between the temperature and the equilibrium 
temperature. The same process is also applied to the electron fraction if we are solving for it. Then the new fluid 
properties are used in another transport and solve iteration. The results presented here are the result of 20 such 
iterations. 

We expect that the equilibrium temperature Tgq and electron fraction Yg^eq should settle to values such that 
Mr/e,EOs(p,req,y'e,eq) = input and T^q = Tinput, where = l^e + fJ'p - l^n IS given by the EOS. Fig- 

ure 16 demonstrates the estimated and calculated equilibrium temperature and electron fraction at several values of 
1/^5 Tinput 5 Te,inputwhcrc Yg input ^ proxy for the chemical potential inputs, such that /izyg,EOs(P 5 Tinput 5 Te,input) 
input- The equilibrium values of T and Yg are each determined with independent iterative calculations. There 
are regions at low temperatures and high densities where the correct solution is not found, as in the previous test. 
This is a combination of inadequate energy resolution and range in the NuLib tables, as the peaks of the emissivity 
spectra approach either end of the energy range or become too sharp to resolve. When equilibrium temperature and 
electron fraction are calculated simultaneously (not plotted), the now two-dimensional solver is less robust and only 
consistently reaches the correct solution when p < 3 X 1013 g cin-3 and T > 2 MeV. 

C. NEUTRINO PAIR ANNIHILATION 

In each grid cell, Sedonu records the neutrino distribution function by accumulating neutrino energy density in 
bins of neutrino species, neutrino energy, and direction. We use this information to calculate neutrino annihilation 
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Figure 16. Blackbody Irradiation Test: The equilibrium fluid temperature (SolvedT) and electron fraction (Solved Ye) as a function 
of fluid density (Input p), neutrino temperature (Input T), and neutrino chemical potential (via the proxy Input Ye)- The three top 
panels show results where only fluid temperature is solved for, and the three bottom panels show results where only fluid electron 
fraction is solved for. The Sedonu-calculated energy densities (black crosses) match the theoretical ones (red lines) in a wide range of 
fluid conditions. The equilibrium solver has difficulty converging on an electron fraction when the chemical potential is large relative 
to the temperature (i.e. high p or low T), since the neutrino Fermi-Dirac distributions become too sharp to be resolved by the NuLib tables. 


rates. The derivation here assumes that the resulting electrons and positrons are extremely relativistic, i.e., that the 
incoming neutrino energies are much larger than the sum of the electron and positron rest masses. This assumption 
is well justified for most of the neutrino energy range we consider. 

The general neutrino annihilation rate is given by Ruffert et al. (1997), equation 1 , representing the rate of energy 
deposition from the annihilation of neutrinos and anti-neutrinos, which we reproduce here for completeness: 
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Barred quantities refer to the anti-neutrino species, Ej^ is the neutrino energy, rUe is the electron mass, c is the speed 
of light, 0 is the angle between the two incoming neutrinos, and ao = 1.76 x 10 “"^^ cm^ is the fiducial weak interaction 
cross section. The weak coupling constants depend on the neutrino species that is annihilating. For electron neutrino 
and electron anti-neutrino (other species) annihilation, Ci E C 2 ^ 2.34(0.50) and C 3 ^ 1.06 (—0.16). / is the phase 
space distribution function (values lie between 0 and 1 ) and fv^/(? is the neutrino energy density per unit neutrino 
energy per steradian of direction. The latter quantity only differs by a factor of c from the specific intensity ly used 
in Dessart et al. (2009). For simplicity, we perform a first-order numerical integral, and assume the energies and 
directions are confined to delta functions at the energy/direction bin centers. Applying this assumption, we arrive at 
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where Ey^i and Ey^k are energy bin centers for the neutrino and anti-neutrino species, respectively. Cy^ij and Ey^ki are 
the integrated (i.e. total measured in the transport simulation) energy density in the corresponding energy/direction 
bin. Oji is the angle between the centers of direction bins j and 1. This must then be summed over all three neutrino 
species-anti-species pairs to get the total energy deposition rate. 

Recall that we differ slightly from this formalism in that we subtract the outgoing lepton masses from the incoming 
neutrino energies to de-emphasize low-energy neutrino annihilations that do not fit the assumption of high-energy 
neutrinos used in the original derivation. One can subtract the electron masses from the incoming neutrino energies 
and match the result to Equation 11 to see that the annihilation kernel is 
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Table 3 

Results - Global Quantities 
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